Causalidade e Anova

uma introdução

Monty Hall

Paradoxo de Simpson

Paradoxo de Simpson

Paradoxo de Berkson

O que é causalidade?

A pesquisa em causalidade busca descobrir uma relação funcional entre duas ou mais variáveis. Por exemplo:

\[ y = f(x) = a + bx \]

No mundo real, observamos essa relação com ruídos. O importante é que esse ruído não tenha relação com outras variáveis que explicam o fenômeno. Exemplo:

\[ y = f(x) = a + bx + \varepsilon, \ \ \varepsilon\sim\mathcal N(0,\sigma^2) \]

Três passos da causalidade

Vamos pensar em um problema onde estamos estudando se a grama fica molhada, observando o céu e um regador.

  • Predição: 1) dado que eu observei o regador funcionando, qual a probabilidade da grama estar molhada? 2) dado que eu observei o céu nublado, qual a probabilidade da grama estar molhada?
  • Intervenção: 1) se eu girar o regador, qual a probabilidade da grama se molhar? 2) Se eu fizer a dança da chuva, qual a probabilidade da grama se molhar?
  • Contrafactual: 1) Sabendo que eu observei o regador desligado e o chão seco, qual a probabilidade do chão estar molhado se eu tivesse ligado o regador? 2) Sabendo que eu observei o céu azul e o chão seco, qual a probabilidade do chão estar molhado se eu tivesse feito a dança da chuva?

Judea Pearl

Três estruturas básicas

Mediadores

A -> B -> C

. . .

Confundidores / garfos

A <- B -> C

. . .

Colisores

A -> B <- C

Importância em regressão - mediador

lr X1 X1 X2 X2 X1->X2 Y Y X2->Y 
n <- 1000
x1 <- rnorm(n)
x2 <- 1 + x1 + rnorm(n)
y <- 1 + x2 + rnorm(n)
modelo <- lm(y ~ x2)
coef(modelo)
(Intercept)          x2 
  0.9939858   1.0164631 
modelo <- lm(y ~ x1)
coef(modelo)
(Intercept)          x1 
   2.027168    1.022733 
modelo <- lm(y ~ x1 + x2)
coef(modelo)
(Intercept)          x1          x2 
 1.02429874  0.06172073  0.98769209

Importância em regressão - confundidor

lr X1 X1 X2 X2 X1->X2 Y Y X1->Y X2->Y 
n <- 1000
x1 <- rnorm(n)
x2 <- 1 + 2 * x1 + rnorm(n)
y <- 1 + x1 + x2 + rnorm(n)

modelo <- lm(y ~ x1)
coef(modelo)
(Intercept)          x1 
   2.013948    2.941999 
modelo <- lm(y ~ x1 + x2)
coef(modelo)
(Intercept)          x1          x2 
  1.0607408   1.0741558   0.9601979

Importância em regressão - colisor

lr X2 X2 X1 X1 X2->X1 Y Y Y->X1 
x2 <- rnorm(n)
y <- rnorm(n)
x1 <- 1 + 2 * x2 + 3 * y + rnorm(n)
lm(y ~ x2) |>
  summary() |>
  broom::tidy()
# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)  -0.0148    0.0327    -0.452   0.651
2 x2            0.0350    0.0337     1.04    0.300
lm(y[x1 > 1] ~ x2[x1 > 1]) |>
  summary() |>
  broom::tidy()
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    0.802    0.0365     22.0  9.94e-75
2 x2[x1 > 1]    -0.302    0.0380     -7.96 1.27e-14

Paradoxos

Monty Hall

Colisor!

lr Escolha1 Escolha1 Apres Apres Escolha1->Apres Carro Carro Carro->Apres

Paradoxo de Berkson

Colisor!

lr Doenca1 Doenca1 Hospital Hospital Doenca1->Hospital Doenca2 Doenca2 Doenca2->Hospital

Paradoxo de Simpson

Paradoxo de Simpson

Confundimento!

lr IdadeSex IdadeSex Tratamento Tratamento IdadeSex->Tratamento Melhorar Melhorar IdadeSex->Melhorar Tratamento->Melhorar

Nesse caso, temos de controlar por Idade e sexo

Outro exemplo!

Confundimento!

lr Pressao Pressao Melhorar Melhorar Pressao->Melhorar Tratamento Tratamento Tratamento->Pressao Tratamento->Melhorar

Nesse caso, não é para controlar por pressão

Modelos conhecidos

Anova

Os passos da ANOVA

  • Definir o que você quer fazer
  • Controle e grupos
  • Definir modelo
  • Testar se existe diferença entre os grupos
  • Comparações múltiplas

As funções da ANOVA

  • lm() ajusta um modelo linear.
  • anova() serve para testar modelos encaixados.
  • aov() é como se fosse anova() + lm() juntos.
  • Para testes múltiplos, utilizamos tukeyHSD().

Qual a ideia do ANOVA?

\[ SQT = SQR + SQE \]

  • \(SQT\): Soma de quadrados total \(\sum_{ij}(y_{ij}-\bar y)^2\)

  • \(SQR\): Soma de quadrados dentro de cada grupo \(\sum_{j}\sum_{i}(y_{ij}-\bar y_j)^2\). Quero que isso seja pequeno.

  • \(SQE\): Soma de quadrados entre grupos \(\sum_{j}(\bar y_{j}-\bar y)^2\). Quero que isso seja grande.

A estatística de teste é dada por

\[ \frac{SQE / gl_E}{SQR / gl_R} \]

Comparando aov() e lm()

aov() testa o fator todo de uma vez.

modelo_aov <- aov(comprimento_bico ~ ilha, data = pinguins)
broom::tidy(modelo_aov)
# A tibble: 2 × 6
  term         df sumsq meansq statistic   p.value
  <chr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>
1 ilha          2 1566.  783.       30.9  4.86e-13
2 Residuals   339 8599.   25.4      NA   NA

lm() testa cada nível do fator individualmente

modelo_lm <- lm(comprimento_bico ~ ilha, data = pinguins)
broom::tidy(modelo_lm)
# A tibble: 3 × 5
  term          estimate std.error statistic   p.value
  <chr>            <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)      45.3      0.390    116.   4.68e-275
2 ilhaDream        -1.09     0.597     -1.83 6.88e-  2
3 ilhaTorgersen    -6.31     0.806     -7.83 6.44e- 14

Como se relacionam?

É possível reproduzir o comportamento de aov() com a utilização de lm(), um modelo nulo e anova()

modelo_lm <- lm(comprimento_bico ~ ilha, data = pinguins)
modelo_lm_nulo <- lm(comprimento_bico ~ 1, data = pinguins)

# Anova serve para comparar modelos
comparacao_modelos <- anova(modelo_lm, modelo_lm_nulo)
comparacao_modelos
Analysis of Variance Table

Model 1: comprimento_bico ~ ilha
Model 2: comprimento_bico ~ 1
  Res.Df     RSS Df Sum of Sq      F   Pr(>F)    
1    339  8598.6                                 
2    341 10164.2 -2   -1565.6 30.862 4.86e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(comparacao_modelos)
# A tibble: 2 × 7
  term                    df.residual    rss    df  sumsq statistic   p.value
  <chr>                         <dbl>  <dbl> <dbl>  <dbl>     <dbl>     <dbl>
1 comprimento_bico ~ ilha         339  8599.    NA    NA       NA   NA       
2 comprimento_bico ~ 1            341 10164.    -2 -1566.      30.9  4.86e-13

Outros comandos

Outros comandos úteis são:

Função Descrição
confint() Intervalo de confiança para os parâmetros
resid() Resíduos do modelo
fitted() Valores ajustados
AIC() Critério de informação de Akaike
model.matrix() Matriz de planejamento (matriz X) do modelo
linearHypotesis() Teste de combinações lineares de parâmetros
vcov() Matriz de variância-covariância dos parâmetros

Fórmulas

Fórmulas

Objetos de classe formula possuem sintaxe muito conveniente para especificar o modelo estatístico que desejamos ajustar. O símbolo que define esses objetos é o ~.

Estrutura:

ajuste <- lm(resposta ~ explicativas)

Então se o objetivo fosse ajustar o modelo

\[ Y_i = \beta_0 + \beta_1X_i + \varepsilon_i, \]

passaríamos ao R a seguinte fórmula

ajuste <- lm(Y ~ X)

Para incorporar mais variáveis usamos o símbolo +. O modelo

\[ Y_i = \beta_0 + \beta_1X_i + \beta_2Z_i + \varepsilon_i, \]

ficaria traduzido como

ajuste <- lm(Y ~ X + Z)

Os operadores * e :

Utilizamos o símbolo * para introduzir os componentes de interação, além dos componentes aditivos.

ajuste <- lm(Y ~ X * Z)

Teoricamente teríamos, para Z contínua, o modelo de regressão

\[ Y_i = \beta_0 + \beta_1X_i + \beta_2Z_i + \beta_3X_i*Z_i + \varepsilon_i, \]

ou, para Z categórica, o modelo de ANCOVA

\[ Y_{ij} = \alpha_j + \beta_jX_{ij} + \varepsilon_{ij}, \]

O operador : faz com que apenas o componente de interação seja incluído no modelo. O modelo

ajuste <- lm(Y ~ X * Z)

é o mesmo que

ajuste <- lm(Y ~ X + Z + X:Z)

Operadores aritméticos

Os operadores aritméticos exercem funções diferentes em fórmulas.

O sinal de + no exemplo induziu em um modelo aditivo em vez de somar X com Z. Para fazer com que eles assumam seus significados aritméticos temos que utilizar a função I().

Exemplo:

ajuste <- lm(Y ~ I(X + Z))

Agora sim o componente I(X + Z) representa a soma de X com Z. Outros exemplos: I(X^2), I(log(X + 1)), I(sqrt(X+Z*5)).

Tabela de símbolos

Símbolo Descrição
+ X inclui a variável X
- X retira a variável X
X * Z inclui X, Z e a interação entre elas
X : Z inclui apenas o componente de interação entre X e Z
(X + Z + W)^2 inclui X, Z, W e as interações 2 a 2
I(X + Z) Função identidade. Inclui uma variável construída pela soma de X com Z
X - 1 Remove o intercepto (regressão passando pela origem)
. O ponto representa ‘todas as demais variáveis’

linearHypothesis e contrastes

Frequentemente temos interesse em saber se parâmetros são diferentes de zero ou se são diferentes entre si. Para isto, costumamos efetuar testes do tipo Wald para combinações lineares dos parâmetros.

Para este fim, a função linearHypothesis() do pacote car faz o trabalho.

Comparações múltiplas

Outra alternativa é utilizar os testes pareados de Tukey. Nesse caso os valores-p já são ajustados.

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = comprimento_bico ~ ilha, data = pinguins)

$ilha
                      diff       lwr        upr     p adj
Dream-Biscoe     -1.089743 -2.495170  0.3156833 0.1628734
Torgersen-Biscoe -6.306505 -8.203279 -4.4097305 0.0000000
Torgersen-Dream  -5.216762 -7.188974 -3.2445487 0.0000000

Tema do Insper para Quarto. Código no GitHub.